package org.cosmocalc;

public class Simpson {
	public static double qsimp(MathFunction func, double a, double b) {
		double ost = -1.0E30;
		double os = -1E30;
		for (int j = 0; j < JMAX; j++) {
			double st = Trapezoid.trapzd(func, a, b, j + 1);
			s = (4.0 * st - ost) / 3.0; // See NumRec eq. 4.2.4
			if (j > 4) // Avoid spurious early convergence
				if (Math.abs(s - os) < EPSILON * Math.abs(os)
						|| (s == 0.0 && os == 0.0)) {
					return s;
				}
			os = s;
			ost = st;
		}
		System.out.println("Too many steps in qsimp");
		return ERR_VAL;
	}

	private static double s;

	public static final double EPSILON = 1.0E-6;

	public static final int JMAX = 50;

	public static final double ERR_VAL = -1E10;

}